home *** CD-ROM | disk | FTP | other *** search
/ 3D GFX / 3D GFX.iso / amiutils / i_l / irit5 / cagd_lib / cbsp_aux.c < prev    next >
C/C++ Source or Header  |  1995-12-30  |  34KB  |  866 lines
  1. /******************************************************************************
  2. * CBsp-Aux.c - Bspline curve auxilary routines.                      *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. /*****************************************************************************
  13. * DESCRIPTION:                                                               M
  14. * Given a Bspline curve - subdivides it into two sub-curves at the given     M
  15. * parametric value.                                                          M
  16. *   Returns pointer to first curve in a list of two subdivided curves.       M
  17. *   The subdivision is achieved by inserting (order-1) knot at the given     M
  18. * parameter value t and spliting the control polygon and knot vector at that M
  19. * location.                                                                  M
  20. *                                                                            *
  21. * PARAMETERS:                                                                M
  22. *   Crv:       To subdivide at parametr value t.                             M
  23. *   t:         Parameter value to subdivide Crv at.                          M
  24. *                                                                            *
  25. * RETURN VALUE:                                                              M
  26. *   CagdCrvStruct *:  A list of the two subdivided curves.                   M
  27. *                                                                            *
  28. * KEYWORDS:                                                                  M
  29. *   BspCrvSubdivAtParam, subdivision, refinement                             M
  30. *****************************************************************************/
  31. CagdCrvStruct *BspCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t)
  32. {
  33.     CagdBType
  34.     NewCrv = FALSE,
  35.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  36.     int i, j, Len, KVLen, Index1, Index2, Mult,
  37.     k = Crv -> Order,
  38.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  39.     CagdRType TMin, TMax, *LOnePts, *ROnePts, *OnePts, *NewKV,
  40.     **Pts, **LPts, **RPts;
  41.     CagdCrvStruct *LCrv, *RCrv;
  42.     BspKnotAlphaCoeffType *A;
  43.  
  44.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  45.     NewCrv = TRUE;
  46.     Crv = CnvrtPeriodic2FloatCrv(Crv);
  47.     }
  48.  
  49.     Len = Crv -> Length;
  50.     KVLen = k + Len;
  51.     Index1 = BspKnotLastIndexL(Crv -> KnotVector, KVLen, t);
  52.     if (Index1 + 1 < k)
  53.     Index1 = k - 1;
  54.     Index2 = BspKnotFirstIndexG(Crv -> KnotVector, KVLen, t);
  55.     if (Index2 > Len)
  56.     Index2 = Len;
  57.     Mult = k - 1 - (Index2 - Index1 - 1);
  58.  
  59.     CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, Crv, TMin, TMax);
  60.  
  61.     LCrv = BspCrvNew(Index1 + 1, k, Crv -> PType);
  62.     RCrv = BspCrvNew(Len - Index2 + k, k, Crv -> PType);
  63.  
  64.     /* Update the new knot vectors. */
  65.     CAGD_GEN_COPY(LCrv -> KnotVector,
  66.          Crv -> KnotVector,
  67.          sizeof(CagdRType) * (Index1 + 1));
  68.     /* Close the knot vector with multiplicity Order: */
  69.     for (j = Index1 + 1; j <= Index1 + k; j++)
  70.     LCrv -> KnotVector[j] = t;
  71.     CAGD_GEN_COPY(&RCrv -> KnotVector[k],
  72.           &Crv -> KnotVector[Index2],
  73.           sizeof(CagdRType) * (Len + k - Index2));
  74.     /* Make sure knot vector starts with multiplicity Order: */
  75.     for (j = 0; j < k; j++)
  76.     RCrv -> KnotVector[j] = t;
  77.  
  78.     /* Now handle the control polygon refinement. */
  79.     Pts = Crv -> Points,
  80.     LPts = LCrv -> Points,
  81.     RPts = RCrv -> Points;
  82.  
  83.     if (Mult > 0) {
  84.     NewKV = IritMalloc(sizeof(CagdRType) * Mult);
  85.     for (i = 0; i < Mult; i++)
  86.         NewKV[i] = (t == TMax ? t - CAGD_DOMAIN_EPSILON : t);
  87.     A = BspKnotEvalAlphaCoefMerge(k, Crv -> KnotVector, Len, NewKV,
  88.                       Mult);
  89.     IritFree((VoidPtr) NewKV);
  90.     }
  91.     else {
  92.     A = BspKnotEvalAlphaCoef(k, Crv -> KnotVector, Len,
  93.                     Crv -> KnotVector, Len);
  94.     }
  95.  
  96.     /* Note that Mult can be negative in cases where original       */
  97.     /* multiplicity was order or more and we need to compensate     */
  98.     /* here, since Alpha matrix will be just a unit matrix then.    */
  99.     Mult = Mult >= 0 ? 0 : -Mult;
  100.  
  101.     /* Blend Crv into LCrv. */
  102.     for (j = IsNotRational; j <= MaxCoord; j++) {
  103.     LOnePts = LPts[j];
  104.     OnePts = Pts[j];
  105.     for (i = 0; i < LCrv -> Length; i++, LOnePts++)
  106.         CAGD_ALPHA_BLEND(A, i, OnePts, Crv -> Length, LOnePts);
  107.     }
  108.  
  109.     /* Blend Crv into RCrv. */
  110.     for (j = IsNotRational; j <= MaxCoord; j++) {
  111.     ROnePts = RPts[j];
  112.     OnePts = Pts[j];
  113.     for (i = LCrv -> Length - 1 + Mult;
  114.          i < LCrv -> Length + RCrv -> Length - 1 + Mult;
  115.          i++, ROnePts++)
  116.         CAGD_ALPHA_BLEND(A, i, OnePts, Crv -> Length, ROnePts);
  117.     }
  118.  
  119.     BspKnotFreeAlphaCoef(A);
  120.  
  121.     BspKnotMakeRobustKV(RCrv -> KnotVector,
  122.             RCrv -> Order + RCrv -> Length);
  123.     BspKnotMakeRobustKV(LCrv -> KnotVector,
  124.             LCrv -> Order + LCrv -> Length);
  125.  
  126.     LCrv -> Pnext = RCrv;
  127.  
  128.     if (NewCrv)
  129.     CagdCrvFree(Crv);
  130.  
  131.     return LCrv;
  132. }
  133.  
  134. /*****************************************************************************
  135. * DESCRIPTION:                                                               M
  136. *  Inserts n knot, all with the value t. In no case will the multiplicity of M
  137. * a knot be greater or equal to the curve order.                 M
  138. *                                                                            *
  139. * PARAMETERS:                                                                M
  140. *   Crv:        To refine by insertion (upto) n knot of value t.             M
  141. *   t:          Parameter value of new knot to insert.                       M
  142. *   n:          Maximum number of times t should be inserted.                M
  143. *                                                                            *
  144. * RETURN VALUE:                                                              M
  145. *   CagdCrvStruct *:   Refined Crv with n knots of value t.                  M
  146. *                                                                            *
  147. * KEYWORDS:                                                                  M
  148. *   BspCrvKnotInsertNSame, refinement, subdivision                           M
  149. *****************************************************************************/
  150. CagdCrvStruct *BspCrvKnotInsertNSame(CagdCrvStruct *Crv, CagdRType t, int n)
  151. {
  152.     int i,
  153.     CrntMult = BspKnotFindMult(Crv -> KnotVector, Crv -> Order,
  154.                                Crv -> Length, t),
  155.     Mult = MIN(n, Crv -> Order - CrntMult - 1);
  156.     CagdCrvStruct *RefinedCrv;
  157.  
  158.     if (Mult > 0) {
  159.     CagdRType
  160.         *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  161.  
  162.     for (i = 0; i < Mult; i++)
  163.         NewKV[i] = t;
  164.  
  165.     RefinedCrv = BspCrvKnotInsertNDiff(Crv, FALSE, NewKV, Mult);
  166.  
  167.     IritFree((VoidPtr) NewKV);
  168.     }
  169.     else {
  170.     RefinedCrv = CagdCrvCopy(Crv);
  171.     }
  172.  
  173.     return RefinedCrv;
  174. }
  175.  
  176. /*****************************************************************************
  177. * DESCRIPTION:                                                               M
  178. *  Inserts n knot with different values as defined by the vector t. If,      M
  179. * however, Replace is TRUE, the knot are simply replacing the current knot   M
  180. * vector.                                     M
  181. *                                                                            *
  182. * PARAMETERS:                                                                M
  183. *   Crv:        To refine by insertion (upto) n knot of value t.             M
  184. *   Replace:    if TRUE, the n knots in t should replace the knot vector     M
  185. *               of size n of Crv. Sizes must match. If False, n new knots    M
  186. *               as defined by t will be introduced into Crv.                 M
  187. *   t:          New knots to introduce/replace knot vector of Crv.           M
  188. *   n:          Size of t.                                                   M
  189. *                                                                            *
  190. * RETURN VALUE:                                                              M
  191. *   CagdCrvStruct *:   Refined Crv with n new knots.                         M
  192. *                                                                            *
  193. * KEYWORDS:                                                                  M
  194. *   BspCrvKnotInsertNDiff, refinement, subdivision                           M
  195. *****************************************************************************/
  196. CagdCrvStruct *BspCrvKnotInsertNDiff(CagdCrvStruct *Crv,
  197.                      CagdBType Replace,
  198.                      CagdRType *t, int n)
  199. {
  200.     CagdBType
  201.     NewCrv = FALSE,
  202.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  203.     CagdRType *KnotVector;
  204.     int Length,
  205.     Order = Crv -> Order,
  206.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  207.     CagdCrvStruct *RefinedCrv;
  208.  
  209.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  210.     NewCrv = TRUE;
  211.     Crv = CnvrtPeriodic2FloatCrv(Crv);
  212.     }
  213.  
  214.     Length = Crv -> Length;
  215.     KnotVector = Crv -> KnotVector;
  216.  
  217.     if (Replace) {
  218.     int i;
  219.  
  220.     if (Order + Length != n)
  221.         CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
  222.     for (i = 1; i < n; i++)
  223.         if (t[i] < t[i - 1])
  224.         CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
  225.  
  226.         RefinedCrv = CagdCrvCopy(Crv);
  227.     for (i = 0; i < n; i++)
  228.         RefinedCrv -> KnotVector[i] = *t++;
  229.     }
  230.     else if (n == 0) {
  231.     RefinedCrv = CagdCrvCopy(Crv);
  232.     }
  233.     else {
  234.     int i, j, LengthKVt;
  235.     BspKnotAlphaCoeffType *A;
  236.     CagdRType *MergedKVt, TMin, TMax;
  237.  
  238.     BspCrvDomain(Crv, &TMin, &TMax);
  239.  
  240.     for (i = 1; i < n; i++)
  241.         if (t[i] < t[i - 1])
  242.         CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
  243.     for (i = 0; i < n; i++)
  244.         if (t[i] >= TMax)
  245.             t[i] -= CAGD_DOMAIN_EPSILON;
  246.  
  247.     /* Compute the Alpha refinement matrix. */
  248.     MergedKVt = BspKnotMergeTwo(KnotVector, Length + Order,
  249.                     t, n, 0, &LengthKVt);
  250.     A = BspKnotEvalAlphaCoef(Order, KnotVector, Length,
  251.                  MergedKVt, LengthKVt - Order);
  252.  
  253.     RefinedCrv = BspCrvNew(Crv -> Length + n, Order, Crv -> PType);
  254.  
  255.     /* Update the knot vector. */
  256.     IritFree((VoidPtr) RefinedCrv -> KnotVector);
  257.     RefinedCrv -> KnotVector = MergedKVt;
  258.  
  259.     /* Update the control mesh */
  260.     for (j = IsNotRational; j <= MaxCoord; j++) {
  261.         CagdRType
  262.         *ROnePts = RefinedCrv -> Points[j],
  263.         *OnePts = Crv -> Points[j];
  264.         for (i = 0; i < RefinedCrv -> Length; i++, ROnePts++)
  265.         CAGD_ALPHA_BLEND(A, i, OnePts, Crv -> Length, ROnePts);
  266.     }
  267.     BspKnotFreeAlphaCoef(A);
  268.     }
  269.  
  270.     BspKnotMakeRobustKV(RefinedCrv -> KnotVector,
  271.             RefinedCrv -> Order + RefinedCrv -> Length);
  272.  
  273.     if (NewCrv)
  274.     CagdCrvFree(Crv);
  275.  
  276.     return RefinedCrv;
  277. }
  278.  
  279. /*****************************************************************************
  280. * DESCRIPTION:                                                               M
  281. * Returns a new curve, identical to the original but with order N.         M
  282. *   Degree raise is computed by multiplying by a constant 1 curve of order   M
  283. *                                                                            *
  284. * PARAMETERS:                                                                M
  285. *   Crv:        To raise its degree to a NewOrder.                           M
  286. *   NewOrder:   NewOrder for Crv.                                            M
  287. *                                                                            *
  288. * RETURN VALUE:                                                              M
  289. *   CagdCrvStruct *:  A curve of order NewOrder representing the same        M
  290. *                     geometry as Crv.                                       M
  291. *                                                                            *
  292. * KEYWORDS:                                                                  M
  293. *   BspCrvDegreeRaiseN, degree raising                                       M
  294. *****************************************************************************/
  295. CagdCrvStruct *BspCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder)
  296. {
  297.     CagdBType
  298.     NewCrv = FALSE;
  299.     int i, j, RaisedOrder, Length, KvLen1,
  300.     Order = Crv -> Order,
  301.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  302.     CagdCrvStruct *RaisedCrv, *UnitCrv;
  303.     CagdRType *Kv;
  304.  
  305.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  306.     NewCrv = TRUE;
  307.     Crv = CnvrtPeriodic2FloatCrv(Crv);
  308.     }
  309.  
  310.     Length = Crv -> Length;
  311.     KvLen1 = Order + Length - 1;
  312.     Kv = Crv -> KnotVector;
  313.  
  314.     if (NewOrder < Order) {
  315.     CAGD_FATAL_ERROR(CAGD_ERR_WRONG_ORDER);
  316.     return NULL;
  317.     }
  318.     RaisedOrder = NewOrder - Order + 1;
  319.  
  320.     UnitCrv = BspCrvNew(RaisedOrder, RaisedOrder,
  321.             CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  322.     for (i = 0; i < RaisedOrder * 2; i++)
  323.     UnitCrv -> KnotVector[i] = i >= RaisedOrder ? Kv[KvLen1] : Kv[0];
  324.     for (i = 1; i <= MaxCoord; i++)
  325.     for (j = 0; j < RaisedOrder; j++)
  326.         UnitCrv -> Points[i][j] = 1.0;
  327.  
  328.     RaisedCrv = BspCrvMult(Crv, UnitCrv);
  329.  
  330.     CagdCrvFree(UnitCrv);
  331.  
  332.     if (NewCrv)
  333.     CagdCrvFree(Crv);
  334.  
  335.     return RaisedCrv;
  336. }
  337.  
  338. /*****************************************************************************
  339. * DESCRIPTION:                                                               M
  340. * Returns a new curve, identical to the original but with one degree higher. M
  341. *                                                                            *
  342. * PARAMETERS:                                                                M
  343. *   Crv:        To raise it degree by one.                                   M
  344. *                                                                            *
  345. * RETURN VALUE:                                                              M
  346. *   CagdCrvStruct *:  A curve with one degree higher representing the same   M
  347. *                     geometry as Crv.                                       M
  348. *                                                                            *
  349. * KEYWORDS:                                                                  M
  350. *   BspCrvDegreeRaise, degree raising                                        M
  351. *****************************************************************************/
  352. CagdCrvStruct *BspCrvDegreeRaise(CagdCrvStruct *Crv)
  353. {
  354.     CagdBType
  355.     NewCrv = FALSE,
  356.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  357.     int i, i2, j, RaisedLen, Length,
  358.     Order = Crv -> Order,
  359.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  360.     CagdCrvStruct *RaisedCrv;
  361.  
  362.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  363.     NewCrv = TRUE;
  364.     Crv = CnvrtPeriodic2FloatCrv(Crv);
  365.     }
  366.  
  367.     Length = Crv -> Length;
  368.  
  369.     if (Order > 2)
  370.     return BspCrvDegreeRaiseN(Crv, Order + 1);
  371.  
  372.     /* If curve is linear, degree raising means basically to increase the    */
  373.     /* knot multiplicity of each segment by one and add a middle point for   */
  374.     /* each such segment.                             */
  375.     RaisedLen = Length * 2 - 1;
  376.     RaisedCrv = BspCrvNew(RaisedLen, Order + 1, Crv -> PType);
  377.  
  378.     /* Update the control polygon; */
  379.     for (j = IsNotRational; j <= MaxCoord; j++)             /* First point. */
  380.     RaisedCrv -> Points[j][0] = Crv -> Points[j][0];
  381.  
  382.     for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  383.     for (j = IsNotRational; j <= MaxCoord; j++) {
  384.         RaisedCrv -> Points[j][i2] =
  385.         Crv -> Points[j][i-1] * 0.5 + Crv -> Points[j][i] * 0.5;
  386.         RaisedCrv -> Points[j][i2 + 1] = Crv -> Points[j][i];
  387.     }
  388.  
  389.     /* Update the knot vector. */
  390.     for (i = 0; i < 3; i++)
  391.     RaisedCrv -> KnotVector[i] = Crv -> KnotVector[0];
  392.  
  393.     for (i = 2, j = 3; i < Length; i++, j += 2)
  394.     RaisedCrv -> KnotVector[j] = RaisedCrv -> KnotVector[j + 1] = 
  395.         Crv -> KnotVector[i];
  396.  
  397.     for (i = j; i < j + 3; i++)
  398.     RaisedCrv -> KnotVector[i] = Crv -> KnotVector[Length];
  399.  
  400.     if (NewCrv)
  401.     CagdCrvFree(Crv);
  402.  
  403.     return RaisedCrv;
  404. }
  405.  
  406. /*****************************************************************************
  407. * DESCRIPTION:                                                               M
  408. * Returns a unit vector, equal to the tangent to Crv at parameter value t.   M
  409. *   Algorithm: insert (order - 1) knots and return control polygon tangent.  M
  410. *                                                                            *
  411. * PARAMETERS:                                                                M
  412. *   Crv:       Crv for which to compute a unit tangent.                      M
  413. *   t:         The parameter at which to compute the unit tangent.           M
  414. *                                                                            *
  415. * RETURN VALUE:                                                              M
  416. *   CagdVecStruct *:  A pointer to a static vector holding the tangent       M
  417. *                     information.                                           M
  418. *                                                                            *
  419. * KEYWORDS:                                                                  M
  420. *   BspCrvTangent, tangent                                                   M
  421. *****************************************************************************/
  422. CagdVecStruct *BspCrvTangent(CagdCrvStruct *Crv, CagdRType t)
  423. {
  424.     static CagdVecStruct P2;
  425.     CagdVecStruct P1, *T;
  426.     CagdRType TMin, TMax;
  427.     CagdCrvStruct
  428.     *FCrv = CAGD_IS_PERIODIC_CRV(Crv) ? CnvrtPeriodic2FloatCrv(Crv) : Crv;
  429.     int k = FCrv -> Order,
  430.     Len = FCrv -> Length,
  431.     OpenEnd = BspCrvHasOpenEC(FCrv),
  432.     Index = BspKnotLastIndexL(FCrv -> KnotVector, k + Len, t);
  433.     CagdPointType
  434.     PType = FCrv -> PType;
  435.     CagdCrvStruct *RefinedCrv;
  436.  
  437.     CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, FCrv, TMin, TMax);
  438.  
  439.     if (APX_EQ(t, TMin) && OpenEnd) {
  440.     /* Use Crv starting tangent direction. */
  441.     CagdCoerceToE3(P1.Vec, FCrv -> Points, 0, PType);
  442.     CagdCoerceToE3(P2.Vec, FCrv -> Points, 1, PType);
  443.     }
  444.     else if (APX_EQ(t, TMax) && OpenEnd) {
  445.     /* Use Crv ending tangent direction. */
  446.     CagdCoerceToE3(P1.Vec, FCrv -> Points, Len - 2, PType);
  447.     CagdCoerceToE3(P2.Vec, FCrv -> Points, Len - 1, PType);
  448.     }
  449.     else {
  450.     RefinedCrv = BspCrvKnotInsertNSame(FCrv, t, k - 1);
  451.  
  452.     CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
  453.     CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
  454.  
  455.     CagdCrvFree(RefinedCrv);
  456.     }
  457.  
  458.     CAGD_SUB_VECTOR(P2, P1);
  459.  
  460.     if (CAGD_LEN_VECTOR(P2) < IRIT_EPSILON) {
  461.     if (AttrGetIntAttrib(Crv -> Attr, "_tan") != TRUE) {
  462.         /* Try to move a little. This location has zero speed. However,  */
  463.         /* do it only once since we can be here forever. The "_tan"      */
  464.         /* attribute guarantee we will try to move EPSILON only once.    */
  465.         AttrSetIntAttrib(&Crv -> Attr, "_tan", TRUE);
  466.  
  467.         T = BspCrvTangent(Crv, t - TMin < TMax - t ? t + EPSILON
  468.                                : t - EPSILON);
  469.  
  470.         AttrFreeOneAttribute(&Crv -> Attr, "_tan");
  471.         if (FCrv != Crv)
  472.         CagdCrvFree(FCrv);
  473.         return T;
  474.     }
  475.     else {
  476.         /* A zero length vector signals failure to compute tangent. */
  477.         if (FCrv != Crv)
  478.         CagdCrvFree(FCrv);
  479.         return &P2;
  480.     }
  481.     }
  482.     else {
  483.     CAGD_NORMALIZE_VECTOR(P2);            /* Normalize the vector. */
  484.     if (FCrv != Crv)
  485.         CagdCrvFree(FCrv);
  486.     return &P2;
  487.     }
  488. }
  489.  
  490. /*****************************************************************************
  491. * DESCRIPTION:                                                               M
  492. * Returns a unit vector, equal to the binormal to Crv at parameter value t.  M
  493. *   Algorithm: insert (order - 1) knots and using 3 consecutive control      M
  494. * points at the refined location (p1, p2, p3), compute to binormal to be the M
  495. * cross product of the two vectors (p1 - p2) and (p2 - p3).             M
  496. *   Since a curve may have not BiNormal at inflection points or if the 3     M
  497. * points are colinear, NULL will be returned at such cases.             M
  498. *                                                                            *
  499. * PARAMETERS:                                                                M
  500. *   Crv:       Crv for which to compute a unit binormal.                     M
  501. *   t:         The parameter at which to compute the unit binormal.          M
  502. *                                                                            *
  503. * RETURN VALUE:                                                              M
  504. *   CagdVecStruct *:  A pointer to a static vector holding the binormal      M
  505. *                     information.                                           M
  506. *                                                                            *
  507. * KEYWORDS:                                                                  M
  508. *   BspCrvBiNormal, binormal                                                 M
  509. *****************************************************************************/
  510. CagdVecStruct *BspCrvBiNormal(CagdCrvStruct *Crv, CagdRType t)
  511. {
  512.     static CagdVecStruct P3;
  513.     CagdVecStruct P1, P2;
  514.     CagdRType TMin, TMax;
  515.     CagdCrvStruct
  516.     *FCrv = CAGD_IS_PERIODIC_CRV(Crv) ? CnvrtPeriodic2FloatCrv(Crv) : Crv;
  517.     int k = FCrv -> Order,
  518.     Len = FCrv -> Length,
  519.     OpenEnd = BspCrvHasOpenEC(FCrv),
  520.     Index = BspKnotLastIndexL(FCrv -> KnotVector, k + Len, t);
  521.     CagdPointType
  522.     PType = FCrv -> PType;
  523.     CagdCrvStruct *RefinedCrv;
  524.  
  525.     CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, FCrv, TMin, TMax);
  526.  
  527.     /* Can not compute for linear curves. */
  528.     if (k <= 2)
  529.     return NULL;
  530.  
  531.     /* For planar curves, B is trivially the Z axis. */
  532.     if (CAGD_NUM_OF_PT_COORD(FCrv -> PType) == 2) {
  533.     P3.Vec[0] = P3.Vec[1] = 0.0;
  534.     P3.Vec[2] = 1.0;
  535.     return &P3;
  536.     }
  537.  
  538.     if (APX_EQ(t, TMin) && OpenEnd) {
  539.     /* Use Crv starting tangent direction. */
  540.     CagdCoerceToE3(P1.Vec, FCrv -> Points, 0, PType);
  541.     CagdCoerceToE3(P2.Vec, FCrv -> Points, 1, PType);
  542.     CagdCoerceToE3(P3.Vec, FCrv -> Points, 2, PType);
  543.     }
  544.     else if (APX_EQ(t, TMax) && OpenEnd) {
  545.     /* Use Crv ending tangent direction. */
  546.     CagdCoerceToE3(P1.Vec, FCrv -> Points, Len - 3, PType);
  547.     CagdCoerceToE3(P2.Vec, FCrv -> Points, Len - 2, PType);
  548.     CagdCoerceToE3(P3.Vec, FCrv -> Points, Len - 1, PType);
  549.     }
  550.     else {
  551.     RefinedCrv = BspCrvKnotInsertNSame(FCrv, t, k - 1);
  552.  
  553.     CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
  554.     CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
  555.     CagdCoerceToE3(P3.Vec, RefinedCrv -> Points, Index + 2, PType);
  556.  
  557.     CagdCrvFree(RefinedCrv);
  558.     }
  559.  
  560.     CAGD_SUB_VECTOR(P1, P2);
  561.     CAGD_SUB_VECTOR(P2, P3);
  562.  
  563.     CROSS_PROD(P3.Vec, P1.Vec, P2.Vec);
  564.  
  565.     if (FCrv != Crv)
  566.     CagdCrvFree(FCrv);
  567.  
  568.     if ((t = CAGD_LEN_VECTOR(P3)) < IRIT_EPSILON)
  569.     return NULL;
  570.     else
  571.     CAGD_DIV_VECTOR(P3, t);                /* Normalize the vector. */
  572.  
  573.     return &P3;
  574. }
  575.  
  576. /*****************************************************************************
  577. * DESCRIPTION:                                                               M
  578. * Returns a unit vector, equal to the normal of Crv at parameter value t.    M
  579. *   Algorithm: returns the cross product of the curve tangent and binormal.  M
  580. *                                                                            *
  581. * PARAMETERS:                                                                M
  582. *   Crv:       Crv for which to compute a unit normal.                       M
  583. *   t:         The parameter at which to compute the unit normal.            M
  584. *                                                                            *
  585. * RETURN VALUE:                                                              M
  586. *   CagdVecStruct *:  A pointer to a static vector holding the normal        M
  587. *                     information.                                           M
  588. *                                                                            *
  589. * KEYWORDS:                                                                  M
  590. *   BspCrvNoraml, normal                                                     M
  591. *****************************************************************************/
  592. CagdVecStruct *BspCrvNormal(CagdCrvStruct *Crv, CagdRType t)
  593. {
  594.     static CagdVecStruct N, *T, *B;
  595.  
  596.     T = BspCrvTangent(Crv, t);
  597.     B = BspCrvBiNormal(Crv, t);
  598.  
  599.     if (T == NULL || B == NULL)
  600.     return NULL;
  601.  
  602.     CROSS_PROD(N.Vec, T -> Vec, B -> Vec);
  603.  
  604.     CAGD_NORMALIZE_VECTOR(N);                /* Normalize the vector. */
  605.  
  606.     return &N;
  607. }
  608.  
  609. /*****************************************************************************
  610. * DESCRIPTION:                                                               M
  611. * Returns a new curve, equal to the given curve, differentiated once.        M
  612. *   Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M
  613. * Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2.  V
  614. *                                                                            *
  615. * PARAMETERS:                                                                M
  616. *   Crv:        To differentiate.                                            M
  617. *                                                                            *
  618. * RETURN VALUE:                                                              M
  619. *   CagdCrvStruct *:   Differentiated curve.                                 M
  620. *                                                                            *
  621. * KEYWORDS:                                                                  M
  622. *   BspCrvDerive, derivatives                                                M
  623. *****************************************************************************/
  624. CagdCrvStruct *BspCrvDerive(CagdCrvStruct *Crv)
  625. {
  626.     CagdBType
  627.     NewCrv = FALSE,
  628.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  629.     int i, j, Len,
  630.     k = Crv -> Order,
  631.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  632.     CagdRType *Kv;
  633.     CagdCrvStruct *DerivedCrv;
  634.  
  635.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  636.     NewCrv = TRUE;
  637.     Crv = CnvrtPeriodic2FloatCrv(Crv);
  638.     }
  639.  
  640.     if (!IsNotRational) {
  641.     DerivedCrv = BspCrvDeriveRational(Crv);
  642.     if (NewCrv)
  643.         CagdCrvFree(Crv);
  644.     return DerivedCrv;
  645.     }
  646.  
  647.     Len = Crv -> Length;
  648.     Kv = Crv -> KnotVector;
  649.     DerivedCrv = BspCrvNew(MAX(Len - 1, 1), MAX(k - 1, 1), Crv -> PType);
  650.  
  651.     if (k >= 2) {
  652.     for (i = 0; i < Len - 1; i++) {
  653.         CagdRType
  654.         Denom = Kv[i + k] - Kv[i + 1];
  655.  
  656.         for (j = IsNotRational; j <= MaxCoord; j++)
  657.         DerivedCrv -> Points[j][i] = (k - 1) *
  658.             (Crv -> Points[j][i + 1] - Crv -> Points[j][i]) / Denom;
  659.     }
  660.     }
  661.     else {
  662.     for (i = 0; i < MAX(Len - 1, 1); i++)
  663.         for (j = IsNotRational; j <= MaxCoord; j++)
  664.         DerivedCrv -> Points[j][i] = 0.0;
  665.          
  666.     }
  667.  
  668.     CAGD_GEN_COPY(DerivedCrv -> KnotVector,
  669.           &Crv -> KnotVector[k < 2 ? 0 : 1],
  670.           sizeof(CagdRType) * (MAX(k - 1, 1) + MAX(Len - 1, 1)));
  671.  
  672.     if (NewCrv)
  673.     CagdCrvFree(Crv);
  674.  
  675.     return DerivedCrv;
  676. }
  677.  
  678. /*****************************************************************************
  679. * DESCRIPTION:                                                               M
  680. * Returns a new Bspline curve, equal to the integral of the given Bspline    M
  681. * crv.                                                                     M
  682. * The given Bspline curve should be nonrational.                 M
  683. *                                         V
  684. *           l           l           l       l+1                 V
  685. *   /         /-           -      /       -   P    -                 V
  686. *  |        | \        n       \     |  n       \    i   \              n+1         V
  687. *  | C(t) = | / P  B (t) = / P   | B (t) = / -----  / ( t   - t ) B   (t) =  V
  688. * /        /  -     i  i       -  i /   i       - n + 1  -    j+n   j   j         V
  689. *            i=0          i=0             i=0     j=i+1                 V
  690. *                                         V
  691. *        l+1 j-1                                 V
  692. *         -   -                                     V
  693. *     1   \   \            n+1                         V
  694. * = ----- /   / P  ( t   - t ) B   (t)                         V
  695. *   n + 1 -   -  i    j+n   j    j                         V
  696. *        j=1 i=0                                 V
  697. *                                         M
  698. *                                                                            *
  699. * PARAMETERS:                                                                M
  700. *   Crv:         Curve to integrate.                                         M
  701. *                                                                            *
  702. * RETURN VALUE:                                                              M
  703. *   CagdCrvStruct *:   Integrated curve.                                     M
  704. *                                                                            *
  705. * KEYWORDS:                                                                  M
  706. *   BspCrvIntegrate, integrals                                               M
  707. *****************************************************************************/
  708. CagdCrvStruct *BspCrvIntegrate(CagdCrvStruct *Crv)
  709. {
  710.     CagdBType
  711.     NewCrv = FALSE;
  712.     int i, j, k, Len,
  713.     n = Crv -> Order,
  714.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  715.     CagdCrvStruct *IntCrv;
  716.     CagdRType *Kv;
  717.  
  718.     if (CAGD_IS_PERIODIC_CRV(Crv)) {
  719.     NewCrv = TRUE;
  720.     Crv = CnvrtPeriodic2FloatCrv(Crv);
  721.     }
  722.  
  723.     if (CAGD_IS_RATIONAL_CRV(Crv))
  724.     CAGD_FATAL_ERROR(CAGD_ERR_RATIONAL_NO_SUPPORT);
  725.  
  726.     Len = Crv -> Length;
  727.     Kv = Crv -> KnotVector;
  728.  
  729.     IntCrv = BspCrvNew(Len + 1, n + 1, Crv -> PType);
  730.  
  731.     /* Copy the knot vector and duplicate the two end knots. */
  732.     CAGD_GEN_COPY(&IntCrv -> KnotVector[1], Kv, sizeof(CagdRType) * (Len + n));
  733.     IntCrv -> KnotVector[0] = Kv[0];
  734.     IntCrv -> KnotVector[Len + n + 1] = Kv[Len + n - 1];
  735.     Kv = IntCrv -> KnotVector;
  736.  
  737.     for (k = 1; k <= MaxCoord; k++) {
  738.     CagdRType
  739.         *Points = Crv -> Points[k],
  740.         *IntPoints = IntCrv -> Points[k];
  741.  
  742.     for (j = 0; j < Len + 1; j++) {
  743.         IntPoints[j] = 0.0;
  744.         for (i = 0; i < j; i++)
  745.             IntPoints[j] += Points[i] * (Kv[i + n + 1] - Kv[i + 1]);
  746.         IntPoints[j] /= n;
  747.     }
  748.     }
  749.  
  750.     if (NewCrv)
  751.     CagdCrvFree(Crv);
  752.  
  753.     return IntCrv;
  754. }
  755.  
  756. /*****************************************************************************
  757. * DESCRIPTION:                                                               M
  758. * Returns a new linear Bspline curve constructed from the given polyline.    M
  759. *                                                                            *
  760. * PARAMETERS:                                                                M
  761. *   Poly:     To convert to a linear bspline curve.                          M
  762. *                                                                            *
  763. * RETURN VALUE:                                                              M
  764. *   CagdCrvStruct *:  A linear Bspline curve representing Poly.              M
  765. *                                                                            *
  766. * KEYWORDS:                                                                  M
  767. *   CnvrtPolyline2LinBsplineCrv, linear curves, conversion                   M
  768. *****************************************************************************/
  769. CagdCrvStruct *CnvrtPolyline2LinBsplineCrv(CagdPolylineStruct *Poly)
  770. {
  771.     int i,
  772.     Length = Poly -> Length;
  773.     CagdCrvStruct
  774.     *Crv = BspCrvNew(Length, 2, CAGD_PT_E3_TYPE);
  775.     CagdRType
  776.     **Points = Crv -> Points;
  777.     CagdPolylnStruct
  778.     *Pts = Poly -> Polyline;
  779.  
  780.     BspKnotUniformOpen(Length, 2, Crv -> KnotVector);
  781.  
  782.     for (i = 0; i < Length; i++, Pts++) {
  783.     Points[1][i] = Pts -> Pt[0];
  784.     Points[2][i] = Pts -> Pt[1];
  785.     Points[3][i] = Pts -> Pt[2];
  786.     }
  787.  
  788.     return Crv;
  789. }
  790.  
  791. /*****************************************************************************
  792. * DESCRIPTION:                                                               M
  793. * Converts a Bspline curve to a Bspline curve with floating end conditions.  M
  794. *                                                                            *
  795. * PARAMETERS:                                                                M
  796. *   Crv:       Bspline curve to convert to floating end conditions. Assume   M
  797. *              Crv is either periodic or has floating end condition.         M
  798. *                                                                            *
  799. * RETURN VALUE:                                                              M
  800. *   CagdCrvStruct *:  A Bspline curve with floating end conditions,          M
  801. *                     representing the same geometry as Crv.                 M
  802. *                                                                            *
  803. * KEYWORDS:                                                                  M
  804. *   CnvrtPeriodic2FloatCrv, conversion                                       M
  805. *****************************************************************************/
  806. CagdCrvStruct *CnvrtPeriodic2FloatCrv(CagdCrvStruct *Crv)
  807. {
  808.     int i,
  809.     Order = Crv -> Order,
  810.     Length = Crv -> Length,
  811.     MaxAxis = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  812.     CagdCrvStruct *NewCrv;
  813.  
  814.     if (!CAGD_IS_PERIODIC_CRV(Crv)) {
  815.     CAGD_FATAL_ERROR(CAGD_ERR_PERIODIC_EXPECTED);
  816.     return NULL;
  817.     }
  818.  
  819.     NewCrv = BspCrvNew(Length + Order - 1, Order, Crv -> PType);
  820.  
  821.     NewCrv -> KnotVector = BspKnotCopy(Crv -> KnotVector,
  822.                        Length + Order + Order - 1);
  823.  
  824.     for (i = !CAGD_IS_RATIONAL_PT(Crv -> PType); i <= MaxAxis; i++) {
  825.     NewCrv -> Points[i] = (CagdRType *) IritMalloc(sizeof(CagdRType) *
  826.                            (Length + Order + Order - 1));
  827.     CAGD_GEN_COPY(NewCrv -> Points[i], Crv -> Points[i],
  828.               sizeof(CagdRType) * Length);
  829.     CAGD_GEN_COPY(&NewCrv -> Points[i][Length], Crv -> Points[i],
  830.               sizeof(CagdRType) * (Order - 1));
  831.     }
  832.  
  833.     for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
  834.     NewCrv -> Points[i] = NULL;
  835.     
  836.     return NewCrv;
  837. }
  838.  
  839. /*****************************************************************************
  840. * DESCRIPTION:                                                               M
  841. * Converts a Bspline curve to a Bspline curve with open end conditions.      M
  842. *                                                                            *
  843. * PARAMETERS:                                                                M
  844. *   Crv:       Bspline curve to convert to open end conditions.             M
  845. *                                                                            *
  846. * RETURN VALUE:                                                              M
  847. *   CagdCrvStruct *:  A Bspline curve with open end conditions,                 M
  848. *                     representing the same geometry as Crv.                 M
  849. *                                                                            *
  850. * KEYWORDS:                                                                  M
  851. *   CnvrtPeriodic2FloatCrv, conversion                                       M
  852. *****************************************************************************/
  853. CagdCrvStruct *CnvrtFloat2OpenCrv(CagdCrvStruct *Crv)
  854. {
  855.     CagdRType TMin, TMax;
  856.  
  857.     if (!CAGD_IS_BSPLINE_CRV(Crv)) {
  858.     CAGD_FATAL_ERROR(CAGD_ERR_BSP_CRV_EXPECT);
  859.     return NULL;
  860.     }
  861.  
  862.     CagdCrvDomain(Crv, &TMin, &TMax);
  863.     return CagdCrvRegionFromCrv(Crv, TMin, TMax);
  864. }
  865.  
  866.